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ABSTRACT 

We compare the predictions of four different algorithms for the distribution of ionized gas during 
the Epoch of Reionization. These algorithms are all used to run a 100 Mpc/h simulation of reioniza- 
tion with the same initial conditions. Two of the algorithms are state-of-the-art ray-tracing radiative 
transfer codes that use disparate methods to calculate the ionization history. The other two algo- 
rithms are fast but more approximate schemes based on iterative application of a smoothing filter to 
the underlying source and density fields. We compare these algorithms' resulting ionization and 21 cm 
fields using several different statistical measures. The two radiative transfer schemes are in excellent 
agreement with each other (with the cross-correlation coefficient of the ionization fields > 0.8 for 
k < 10 h/Mpc) and in good agreement with the analytic schemes (> 0.6 for k < 1 h/Mpc). When 
used to predict the 21cm power spectrum at different times during reionization, all ionization algo- 
rithms agree with one another at the 10s of percent level. This agreement suggests that the different 
approximations involved in the ray tracing algorithms are sensible and that semi-numerical schemes 
provide a numerically-inexpensive, yet fairly accurate, description of the reionization process. 
Subject headings: cosmology: theory - intergalactic medium - large scale structure of universe 



1. INTRODUCTION 

At z > 20, structures massive enough for the gas to 
cool and form the first galaxies first broke away from the 
Hubble flow. Stars and black holes formed in these struc- 
tures and ultimately ionized and heated the intergalactic 
gas. However, exactly when and how this process oc- 
curred are questions that are still unresolved. There are 
currently two prominent constraints on the reionization 
era. First, Lya forest absorption spectra towards high 
rcdshift quasars show a rapidly increasing opacity to Lya 
photons by the intergalactic medium (IGM) at z > 6 
(e.g. iFan et al.1 120061 ) ). Several studies have interpreted 
this increase as evidence for the end of reionization. At 
the very least, the amount of transmission in the Lyct 
forest at z < 6 indicates that the bulk of reionization 
occurred at higher redshifts. Second, measurements of 
the large-angle cosmic microwave background (CMB) po- 
larization anisotropics suggest that the median redshift 
of the reionization pr ocess was roughly z ~ 10. 4 ± 1.4 
(|Komatsu et al.ll2009fh 

Taken at face value, these obser vations favor an ex- 
tended reionization epoch (see al s o iCenl 120031: IFan et al.1 
[20061 : IBolton fc Haehneltl l2007bl : IWvithe fc Cenl 120071 ). 
However, the interpretation of quasar absorption spectra 
is hampered by the large cross section for Lya absorp- 
tion, which can lead to complete absorption even if the 
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hydrogen is highly ionized. This makes the high-redshift 
Lyq forest spectra difficult to interpret (jLidz et al.ll2006t 
IBecker et al.l 12007] ) . and the data may even be consis- 
tent with reionization completing at redshifts z < 6 
(jLidz et all 120071: iMesingerl l2009j ). Furthermore, the 
CMB polarization measurements onl y offer an integra l 
constraint on th e ionization history (jKogut et al.ll200"3l : 
IPage et al.l[200l . 

Additional model-dependent constraints on reion- 
ization have been derived from other astrophysical 
probes such as (1) the size of th e proximity zone 
around quasars (IWvithe et al.l 120051: IFan et al.l 12006 , 
but see [Mesinger et al.l 120041: IBolton fe Haehneltl 12007a: 
iLidz et all 120071 : iMaselli et al.ll2007t ): (2) a claimed de- 
tection of dampin g wing absorption from neutral IGM 
in quasar spectra (IMesinger fc Haimanl 12004 120071 but 
see lMesinger fe Furlanettoll2008aD : (3) the non-detection 
of intergalactic dam ping wing absorpt i on in a gamma 
ray b urst spectrum (jTotani et al.l 120061 : jMcQuin n et al.l 
l2008f ): and ( 4) the number density and clustering of 
Lyq e mitters (iMalhotra fc Rhoaddl2004t |H aiman fc Ce: 



20051: iFurlanetto et all l2006bt iKashikawa et al 



McQ uinn et al.ll2007aHMesinger fc Furlanettoll2008bl)" 

Redshifted 21 cm emission from the hyperfme tran- 
sition of neutral hydrogen has the potential to pro- 
vide detailed three-dimensional information about the 
evolu tion and morphology o f the reionization process 
(e.g.. Zaldarriasa ct al 2Q0d), Scycr.a! large interfer- 

mmcgjj nirsj i >(.'^ei(\y.i.'c yi, iiuac sg iii.liarva.ijl.edu . ° , 

OiiictCrs win ucVGuu sig,nincanu integration times dur- 
ing the next few years towards detecting this signal. 
These efforts include the M ileura Wide Field Array 
(MWA) (|Bowman et aLll2005[ ) 7 . the Low Frequency Ar- 
ray (LOFAR) 8 t he Giant Metrewave Telescope (GMRT) 
(jPen et al .11 20091) . Precision Array for Probing the Epoch 

7 http:/ /web. haystack. mit.edu/arrays/MWA/ 

8 http://www.lofar.org 
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of Reionization (PAPER) (jParsons et all 12009? ). and 
eventually the Square Kilometer Array (SKA) 9 . 

Our ability to infer information about the prop- 
erties of first sources from the 21cm signal, as well 
as other observations, hinges on the accuracy of 
the theoretical modeling of the ionization struc- 
ture during reionization. Therefore, a number of 
groups have d eveloped 3-D radiative transfer (RT) 
codes (e.g. iGnedinl 120001: iRazoumov et all l200a 
20031: iSokasian et all 120011: iMellema et all 
2007bt ISemelin et all 12 007: 
2008F lAubert fe Tevssierl 



Ciardi et all 
20061: !McC>uinn et all 



Trac fc Cenl 120 07: Alt av et al.1 ._ 

2008t iFinlator et all 120091: iPetkova fc Springe!! 120091 : 
see the recent review in " iTrac fc Cxnedinl 12009ft . Fur- 
thermore, given the large uncertainties in how ionizing 
photons are produced and escape from galaxies (and 
also in the dense systems that by the end of reion- 
ization are the largest sinks of ionizing photons), a 
large region of parameter space needs to be modelled 
in order to interpret observations. These concerns 
have prompted several groups to come up with more 
appro xi mate, but much faster sc heme s ( e.g. iZahn et al 
2001 IMesinger fc Furlanettol [20071 [Geilfc Wvithe 
2008t lAlvarez et all 120091: iChoudhurv et all 120091: 



Thomas et all T2009) . 



Accurate models of the Epoch of Reionization (EoR) 
must include the evolution of the dark matter, gas, radia- 
tion, galaxies, as well as a plethora of feedback processes. 
Although, the morphology of reionizatio n may be robust 
to man y of these modeling uncertainties (McQuin n et all 
I2007H) . Still, the scope of the parameter space is daunt- 
ing. Simulations must resolve the small mass galaxies 
that are expected dominate the ionizing photon budget 
(probably corresponding to the atomic cooling threshold, 
with halo masses of ~ 10 8 M Q at z ~7-10), as well as 
the photon sinks. On the other hand, reionization simu- 
lations must also be large enough to statistically sample 
the distribution of HII regions, which can span tens of 
comovin g megaparsecs in size towards the end of reion- 
ization (IFurlanetto fc Ohl [2001 IZahn et all l200l 120071: 
IMesinger fc Furlanettol 120071 : IShin et all 120081 ). A few 
groups have r ecently come close to achieving this dy- 
namic range (llliev et all l2006at iMcQuinn et all l2007bl: 
ITrac fc Cenll2007l: IShin et al.ll2008t ITrac et alj|2008t see 
the recent review in ITrac fc Gnedinl I2009D . 10 In ad- 
dition, pro mising analytic alternat ives have been sug- 
gested (e.g. IFurlanetto et all [20041 ). The Monte-Carlo 
implementation of these alternatives in large- volume sim- 
ulations allows for a description of the non-spherical 
HII region morphology that can be compared side- 
by-side with simulations (e .g. IZahn et all 120051 120071: 
IMesinger fc Furlanettoll2007l ). 

In this paper, we investigate the convergence of differ- 
ent algorithms for modeling the EoR. We compare two 
different cosmological radiative transfer codes, which fol- 
low ionizing radiation from sources identified as halos 
in an N-body simulation. Our comparison also includes 

9 http://www.skatelescope.org/ 

10 Note however that Lyman limit systems (LLSs), which can 
dominate th e absorption of ionizi ng photons (see for example the 
appendix of Furlanetto & Oh 2005), are still too small to be re- 
solved by state-of-the-art reionization simulations, and therefore 
must be included via some analytic prescription. 



ionization fields generated with somewhat refined ver- 
sions of the fast, semi-numerical algorithms mentioned 
above. The comparisons are done using the same den- 
sity field, and N-body halo field (where applicable). We 
test the codes on the full reionization problem (at least 
to the extent that it is currently simulated) in order to 
study the numerical convergence of the predicted 21 cm 
signal Ou r approach is different but complimentary to 
llliev et all ([2006b. 2009) who compared several RT codes 
using a number of tests problems. Both approaches are 
necessary for studying numerical convergence. 

The structure of this paper is as follows. In Section [2j 
we present the N-body simulation we employ, as well as 
the two different radiative transfer schemes that are run 
in post-processing on top of the N-body simulation. In 
Section[4j we describe the approximate schemes. Section 
[5] investigates the spatial agreement between the ioniza- 
tion fields of these algorithms. Section [6] compares the 
predictions of these schemes for the 21cm power spec- 
trum. In the Appendix, we present further discussion 
and tests of the semi-numerical schemes. 

The simulations and analytic calculations presented 
here are based on a KCDM cosmology with the following 
parameter values: cr 8 = 0.82, h = 0.7, Q m = 0.28, = 
0.046, n s = 0.96, consistent w ith latest constraints from 
WMAP (IKomatsu et alll2009l) . 

2. TEST PROBLEM 

In this paper, the radiative transfer and semi-numeric 
calculations were performed by po st-processing input 
density fields from lTrac et all ((2008). A high- resolution 
N-body simulation with 3072 3 dark matter particles was 
used to evolve the matter distribution in a periodic box of 
comoving volume (100 /i _1 Mpc) 3 . We chose to work with 
82 snapshots of the matter density field (every 10 million 
years from redshift z = 27toz = 6), gridded onto a 
256 3 Cartesian grid for the reionization algorithms. The 
baryons were assumed to trace the dark matter down to 
the grid cell spacing Ax = 390 /i _1 kpc, which is close 
to the anticipated physical smoothing scale for the inter- 
galactic gas (the Jeans length for 10 4 K gas is equal to 
370 h~ 1 kpc at the mean density at z = 6). 

Radiation sources embedded in the large-scale struc- 
ture of the matter distribution were modeled using halos 
catalogued in the N-body simulation. Dark matter halos 
were located on the fly using a friends-of-friends algo- 
rithm, with a linking length of b — 0.2 times the mean in- 
terparticle spacing. Halos with virial temperatures above 
the atomic cooling limit (T > 10 4 K; or M > 10 8 M Q for 
the relevant redshifts) were located with a minimum of 
~ 40 particles (M = 1.0 x 10 8 M Q ), and half of this col- 
lapsed mass budget is resolved with > 400 particles per 
halo. The first source appeared at z ~ 27 within the 
(100 /i _1 Mpc) 3 comoving volume. By z — 6, there were 
> 6 million sources and, when binned on a 256 3 Carte- 
sian grid, the clustered sources occupied ~ 17% of the 
grid. 

3. RADIATIVE TRANSFER ALGORITHMS 

This section describes the two radiative transfer codes 
that are performed on the test problem described in Sec- 
tion [2] the McQuinn et al. code (Section 13- 1[) and the 
Trac & Cen code (Section |3~2| . 
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3.1. McQuinn et al. 
The algorithm prese nted inlMcQuinn et al.l (|2007bT > is 



an improvement of the iSokasian et al.l ( 20011 ) rav tracing 
code. It is a more approximate ray tracing scheme than 
iTrac fe Cenl ((2001 , but has more modest CPU and mem- 
ory requirements. This code was specifically designed to 
perform ray tracing on the millions of sources required 
to study the Eo R and was the first ray tracing code to 
achieve this feat (McQuinn ct al. 2007b). 11 In addition, 
this code has enabled the first param eter space studies 
of this epoch (jMcQuinn et al.ll2007bD . 

This code sends out 768 rays from every source be- 
tween every time step, and rays travel within that time 
step until their photons are expended (i.e., an infinite 
speed of light). It us es the HEALPix-based a daptive ray- 
splitting scheme of lAbel fc Wandelti (|2002t ). such that 
rays split as they travel to ensure that a minimum num- 
ber of rays reach each cell. For the simulations in this 
paper, each cell receives at least 2.1 rays from each source 
when averaged over HEALPix orientations. The num- 
bers 768 (corresponding to HEALPix level 3) and 2.1 
w ere calibrated with conv ergence tests in the Appendix 
of iMcQuinn et al.l (|2007b|) and were the values that this 
code has been run with for large-scale simulations of cos- 
mological reionization. The order that rays are cast in 
this code is randomized over both the sources and direc- 
tions for each time step. Finally, in the test problem in 
this paper, each time step is taken to be At = 10 Myr, 
which corresponds to the time between stored snapshots 
of the density field. Comparable values f or At have been 
used in p ublished studies with this code ([McQuin n et al.l 
I2007allbfi . 

The McQuinn et al. radiative transfer algorithm for 
hydrogen reionization works as follows. At the begin- 
ning of the time step, the ionization fraction in each cell 
is corrected for recombinations that occur during that 
time step. Next, rays are cast. Each ray carries a set 
number of photons and when a ray reaches a cell that 
contains neutral gas, it deposits all of the photons into 
ionizations. The approximation that all the photons are 
absorbed at the ionization front is motivated by the short 
mean free path of ionizing photons in a neutral medium 
for a physically motivated source spectrum during hy- 
drogen reionization (i.e., a mean free path of 10s of kpc 
is expected, which is much smaller than the cell size). 
When the path to the edge of an ionized region becomes 
comparable to the box size, this algorithm slows signifi- 
cantly because it has to follow each ray to the front edge 
each timestep. The solution to the ionization field that 
this algorithm obtains technically depends on the order 
rays are cast from each cell, but convergence studies have 
demonstrated that in practice a nearly i dentical solution 
is reac hed with different ray orderings (McQuin n et al.l 
2007b). After rays have ionized a cell, it is assumed to 
be in photoionization equilibrium with the photoionizing 
background. 

Finally, as we mentioned previously, this algorithm as- 
sumes an infinite speed of light. This approximation is 
the least valid at the end of reionization when the HII 
regions are the largest. Because of this approximation, 

11 Although, HI cooling mass halos were unresolved in the N- 
body simulation used in this study and were instead included in 
post processing prior to ray tracing. 



the McQuinn et al. algorithm produces several per cent 
larger ionization fractions close to the end of reionization, 
compared with the other radiat i ve tra nsfer algorithm in 
this comparison bv ITrac fc Cenl (pOOl . This issue could 
be alleviated in higher resolution simulations than per- 
formed here where Lyman limit absorption systems are 
resolved, which will limit the mean free path. 

3.2. Trac & Cen 

The ITrac fe Cenl (|2007l ) radiative transfer algorithm 
follows the propagation of ionizing photons using an 
adaptive ray tracing technique, featuring several novel 
approaches that distinguish it from other previous ray 
tracing algorithms. The major characteristics include 
conservation of photons, a time-dependent solution that 
is causal, and computational scaling that is independent 
of the number of sources. Here we provide a basic de- 
scription of the algorithm and highlight differences with 
other approaches. 

As a basis, t he a daptive ray-splitting scheme of 
lAbel fe Wandelti (|2002t) is used to efficiently improve spa- 
tial resolution. For a given source, a preset number of 
base rays are cast and they travel a short distance before 
each splits into 4 daughter rays. Successive generations of 
splitting continue downstream such that the angular res- 
olution continually increases farther away from sources. 
This procedure ensures that a minimum number of rays 
intersect each grid cell element of the density field within 
the ionized region. For this comparison study, the num- 
ber of base rays was set at 192 and the minimum number 
of rays intersecting each grid cell was set at 2. 

One major feature of the algorithm is that new source 
rays are cast for every radiative transfer time step, which 
is set equal to the light-crossing time, defined as the 
length of a grid cell divided by the speed of light. Cor- 
respondingly, all rays are advanced a maximum distance 
(when the medium is optically thin) equal to a grid cell 
length in each time step. This synchronization ensures 
a causal solution, a characteristic usually not shared by 
other ray tracing algorithms where rays from multiple 
sources are treated independently. Generally, this ap- 
proach is costly because the light-crossing time at high 
resolution can be quite small compared to the dura- 
tion of reionization. However, Pawlik fc Schave] (2008) 
have pointed out that when the radiation filling factor 
is close to unity near the end of reionization, the light- 
crossing time stepping can actually be computationally 
faster than the infinite speed of light approach. 

Another major feature of the algorithm is an adaptive 
ray-merging scheme designed to circumvent the 0(N 2 ) 
scaling in the multi-source problem. In the brute-force 
approach, every source must emit roughly as many rays 
as there are density grid cells when the radiation filling 
factor is close to unity. For a 256 3 grid with close to 3 
million source cells, a computationally infea sible number 
of 50 t rillion rays would be required in total. ITrac fc Cenl 
(|2007f ) adaptively restrict the maximum number of rays 
entering and exiting a given cell by merging near-parallel 
rays from multiple sources. For this comparison study, 
approximately 100 rays per cell and as many as 2 billion 
rays to propagate within a time step were required by 
the end of reionization. One disadvantage is that a large 
amount of memory is still required to store the neces- 
sary rays. In the current version of the implementation, 
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each ray requires 44 bytes to store intrinsic variables plus 
another 8 bytes per frequency bin to keep track of the 
frequency and photon count. 

For each radiative transfer time step, the photoioniza- 
tion rates are calculated on the same 256 3 grid using 
the incident radiation flux from the rays. The equations 
for the ionization evolution are then solved using a non- 
equilibrium solver where the time integration involves a 
stable implicit scheme. For this comparison study, only 
photoionization at one frequency is considered. However, 
the radiative transfer algorithm can also handle photo- 
heating and multiple fr equencies in a c osmological hy- 
drodynamic simulation (Trac et al. 2008). 

4. SEMI-NUMERIC MODELING 

Our comparisons also include a semi-numerical 
algorithm for generating the ionization field. These 
techniques are motivated by the intuition that sources 
ionize the regions immediately around them and that 
the clustering of sources drives the structure of reion- 
ization. The radiative transfer algorithms described 
previously can require large computer clusters and can 
take days to process a single redshift output. However, 
for many cosmological reionization applications, it might 
be useful to sacrifice some accuracy for an increase in 
speed and dynamic range. This is potentially quite 
useful for many applications, such as in simulating the 
signal for the upcoming 21cm surveys. The effective 
resolution of all planned instruments, including the 
SKA, is much worse than the resolution of the typical 
RT scheme. Additionally the field of view (FoV) 
of some of these instruments will be enormous (e.g. 
5000x5000 Mpc 2 at z ~ 8 for MWA), far out of reach 
of conventional RT algorithms hoping to include sources 
down to the atomic hydrogen cooling threshold 12 . Also, 
very little is presently known about the properties of 
the underlying sources, translating to an enormous 
parameter space in need of exploration. Thus, a "rel- 
atively" accurate yet fast ionization algorithm can be 
an invaluable tool in reionization studies. Below we 
detail the semi-numerical algorithm we use, capable 
of generating large-scale 3D ionization fi elds 13 in a 
matter of minutes on a single processor (see IZahn et al. 
2005 1 [20071: IMesinger fc Furlanettol[2007t IGeil fe Wvithe 
2008t lAlvarez et all 120091: iChoudhurv et al.l 120091: 
Tho mas et al.N2009l for related/alternate semi- numerical 
ionization algorithms) . 

Our semi-numerical algorithm is based on the 
excursio n-set formalism f o r mod eling reionization devel- 
oped by iFurlanetto et al.l (|2004f) . The foundation of this 
approach is to require that the number of ionizing pho- 
tons inside a region be larger than the number of hydro- 
gen atoms it contains. Then ionized regions are identified 

12 It is not clear how useful simulating the FoV of these instru- 
ments would be, given that foreground subtraction will likely ren- 
der most of the transverse mo des unusable for studying the EoR 
(e.g. Morales & Wvithe 2009). However, in order to get proper 
sampling of the line-of-sight (LOS) power, one still needs many 
EoR realizations, as otherwise the LOS modes will add coherently. 
It would be most desirable to correctly include cosmic variance in 
these realizations, either through a very large box/dynamic range, 
or through many smaller boxes where long wavelengths modes have 
been added-in "by hand". 

13 "Ionization fields" is synonymous with "distribution of HII 
regions" in this paper. 



via an excursion-set approach, starting at large scales 
and progressing to small scales, analogous to the deriva- 
tion of the Press-Schecht er mass functions (|Bond et al.l 
flflflltlLacev fc Colelll993l). 

We extend the IFurlanetto et all (|2004h scheme to 3D 
realizations by requiring that an ioniz ed simulation cell 
at coordinate x meet the cri teria (c.f. IZahn et al.l 120051 : 
IMesinger fc Furlanettoll2007D . 

/coll(x,Z,i?) >C' , (1) 

where ^ is some efficiency parameter and / co n is the frac- 
tion of mass residing in collapsed halos inside a sphere 
of mass M = 4/3ttR 3 p[1 + (<5nl).r], with mean phys- 
ical overdensity (<5nl),r, centered on Eulcrian coordi- 
nate x, at redshift z. Following the excursio n-set ap- 
proach (jBond et all Il99lt lLacev fc CoiellT99"l . we de- 
crease the filter scale R 7 starting from some maximum 
value i?max 14 , and progressing down to the cell size, i? C eii, 
with logarithmic filter spacing of i? n oxt = 0.9i? prC v If at 
any filter scale the criterion in eq. Q] is met, this cell is 
flagged as ionized. Note tha t we only flag the central fil- 
ter cell at x as ionized (as in IZahn et al.ll2005l ) instead of 
"painting" the entire sphere w ith radius R as ionized (ala 
IMesinger fc Furlanettol[2007l) . as the former approach is 
faster and thus more versatile. 

Unlike in IZahn et all (|2005D and 
IMesinger fc Furlanettol (|2007h . we take into account 
partially-ionized cells by setting the cell's ionized 
fraction to C/coii(x, z, R cc ii) at the last filter step for 
those cells which are not fully ionized. More specifically, 
we use the unfiltered value from the cubical simulation 
cell instead of filtering it with a spherical filter on 
the cell size, in order to remove aliasing effects from 
using an analytic filter form. A two-phased medium 
(containing only fully ionized and fully neutral cells) is 
a good approximation for stellar sources and simulation 
resolution we study (see the Appendix). Nevertheless, 
partially ionized cells become much more common as 
the cell size is increased. Our prescription accounts for 
cells which have bee n partially ionized by the sources 
they contain (see also lGeil fc W vithe 20Q8|). 

We present two main variations of this semi-numerical 
ionization algorithm, distinguished by their method of 
computing / co n and choice of filter: 

• Fast Fourier Radiative Transfer (FFRT): This vari- 
ant computes / co n using th e conditional Press- 
Schechter (PS) formalism (jLacev fc Cold 119931 : 
ISomerville fc Kolatt] 119991) directly on the evolved 
(i.e., non-linear) density field. Although the col- 
lapse criterion for the conditional Press-Schechter 
was derived using the linear density field, we find 
that using the evolved density from simulations is 
a better match to the RT ionization fields. This is 
to be expected because the evolution leads to non- 

14 flmai can in principle can be set to correspond to 
the mean free path of ionizing photons (Furlanctto & Oh] 120051 ; 
IChoudhurv et al.ll2009T ). As the mean value and spatial fluctua- 
tions of the mean free path are poorly known at present, and as 
they are not fully included in the numerical simulations, for pur- 
poses of comparison we use -R ma x ~ 89 Mpc, which corresponds 
to the entire simulation box. Note however, that in the regime of 
interest when -R max 3> the typical HII region size, the precise value 
of -Rmax does not affect ionization morphology. 
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linear shifts with respect to the linear field. There- 
fore, this scheme bypasses the need to resolve ha- 
los, allowing for a much larger dynamical range. 
The other modification compared to earlier semi- 
numeric algorithms is the use of a sharp k-space 
filter instead of spherical top-hat 15 . 

• Fast Fourier Radiative Transfer- Sources (FFRT- 
S): This variant computes / co ii by filtering the 
N-body halo field (more precisely the star- 
formation rate field). This scheme therefore 
requires a discrete source field, o btained ei- 
ther through N-body simulations (|Zahn et al.l 
120071) or through another sem i-numerical scheme 
(jMesinger fc Furlanettol 120071) . This algorithm by 
default uses a top-hat filter, in order to avoid "ring- 
ing" artifacts in the ionization field which can re- 
sult from the oscillatory configuration-space filter 
response of the sharp k-space 16 . 

In setting the partial ionizations in the FFRT scheme, 
we also account for Poisson fluctuations in the halo 
number, when the mean collapse fraction obtained 
through conditional Press-Schechter 17 becomes small, 
/ co ii(x,z, i? cc ii)[l + <5 cc ii(x, z,R ce n)]M ce ii < 50M min , 
where M ce u is the mean, total (DM+baryonic) mass cor- 
responding to a cell, <5 ce ii is the fractional overdensity in 
the cell's mass, and our faintest ionizing sources corre- 
spond to halo masses of M min ~ 10 8 M Q in this appli- 
cation. We make the simplifying assumption that the 
sub-grid sources are all of the same mass, M m ; n (likely a 
decent approximation for partially ionized cells, since the 
mass function is fairly steep in this regime) . We sample a 
Poisson distribution with mean / C oiiM ce ii(l + S ce \i)/M m i n 
to obtain the number of sub-grid sources, iV m i n . We 
then compute that cell's collapse fraction, f co \\ sampled = 
N m i n M m i n /M ce ii/(l + <Wi)- Poisson fluctuations are 
found to be important when the cell size increases to 
> 1 Mpc (see Fig. |8]and associated discussion). Note 
however that this feature should be left as an option, de- 
pending on the particular application. Although improv- 
ing the statistics of the ionization fields at a fixed red- 
shift, such stochastic fluctuations in the collapsed mass 

15 Although the spherical top-hat is more physically intuitive 
(ionizing photons in a sphere are compared to neutral hydrogen 
atoms in the same sphere), the sharp k-space filter results in a 
somewhat better correlation of the ionization and density fields in 
the FFRT scheme, which is of importance in predicting the 21cm 
signal (see Fig. [7] and associated discussion). Another advantage 
of this new filter is that C/coll wm t> e to equal the global mean 
ionized fraction, in the limit of large box sizes (see Appendix of 
Zahn ct al. 2007). The top-hat filter performs much worse in this 
respect. Therefore, the redshift evolution of reionization is more 
accurately predicted using a sharp k-space filter than using a top- 
hat filter. 

16 In practice, we do not find evidence for such artifacts for the 
regimes studied here. However, they should become more common 
in the regime where sources are rare. Therefore our fiducial scheme 
operating on discrete sources uses a top-hat filter. 

17 S pecifically we use the "hybrid" form of Barkana & Loeb 
(2004), which includes the additional factor o f the ratio of the 
PS and Sheth-Tormcn (Shcth & Tormcn 1999) collapse fractions 
at mean density. Barkana & Loeb (2004) find that this change in 
the overall normalization yields mass functions which better com- 
pare to N-body simulations. Note that, aside from computing the 
partially ionized cells, the difference between this hybrid form and 
PS can be absorbed in the £ efficiency parameter. 



and ionization fields make it impossible to determinis- 
tically track the redshift evolution of a particular re- 
alization. Halos would effectively "pop in-and-out of" 
existence from on e redshift/time step to the next (e.g. 
ISantos et al1[2009h . 

This semi-numerical scheme can be extended to con- 
tain additional physics such as spatially inhomogeneous 
recombinations (jChoudhurv et al.l 120091). a mass d epen- 
dent ionization effic iency (iFurlanetto et al.l l2006a). and 
radiative feedback (jGeil fc Wvithd 120081 ) Since these 
processes are also poorly understood and the RT schemes 
we compare against do not take them into account, we 
will limit our refinements in this work to the simplest 
case that ignores feedback and inhomogenous recombi- 
nations. 

5. COMPARISON OF IONIZATION FIELDS 

This section compares the ionization field s gener - 
ated from the fou r schemes: IMcQuinn et al.l (|2007bj) . 
iTrac fc Cenl (|2007l ). FFRT, and FFRT-S. Figure □ shows 
ionization maps for three different volume- weighted ion- 
ization fractions, xiy = 0.25,0.51, and 0.72. The maps 
are from the same slice, spanning 100 Mpc/h by 100 
Mpc/h with a depth of 0.4 Mpc/h, through the simula- 
tion box. In all of the simulations, the reionization pro- 
cess appears qualitatively similar. In the early stages, 
HII regions typically contain relatively few galaxies and 
are spatially clustered similar to the sources. The HII 
regions expand and merge with one another as reioniza- 
tion progresses. By the time the process is half com- 
pleted, the ionization structures appear much more con- 
nected, resembling the large-scale-structure of the uni- 
verse. The largest HII regions now contain tens of thou- 
sands of galaxies. In the final stages, the ionization fronts 
quickly expand into the remaining neutral regions until 
there is complete overlap. 

The maps from the two RT simulations show the clos- 
est agreement, with very similar ionization structures at 
all times. The ionization fronts are similarly resolved, 
such that even small neutral regions between converging 
fronts are well preserved. Of the two semi-numeric meth- 
ods, the FFRT-S produces maps that are in closer agree- 
ment with the RT results. Although reproducing the 
HII morphology on moderate to large scales, the semi- 
numeric algorithms tend to generate more connected HII 
regions. Such differences are understandable, as the 
semi-numeric models do not include the physics needed 
to exactly capture the position and structure of the ion- 
ization fronts. In the rest of this section, we compare 
several statistics to quantify the reionization morphol- 
ogy and to check the degree of convergence between the 
different schemes. 

5.1. Bubble size distribution 

Figure [5] shows the ionized bubble size distribution, 
defined as the probability that a randomly chosen ion- 
ized cell from the simulation volume is part of an ion- 
ized region of radius R. An ionized cell is assigned to 
be in a bubble of radius R, where R is set equal to the 
largest sphere surroundin g the cell whose m ean enclosed 
ionized fraction is 0.9 (see lZahn et al.ll2007l for more dis- 
cussion). This definition is convenient because it is sim- 
ilar to t hat which is used in e xcursion set reionization 
models (|Furlanetto et al.l 120041 ) . In addition, it has the 
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z=8.49 z=7.56 z=7.11 

X=0.25 X=0.51 X=0.72 




Fig. 1. — Comparison of ionization fields generated from four schemes: McQuinn et al., Trac & Cen, FFRT-S, and FFRT. The maps are 
from the same slice (100 Mpc/h by 100 Mpc/h with depth of 0.4 Mpc/h) through the simulation box. 
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Fig. 2. — The ionized bubble size distribution for the four 
schemes: Trac & Cen (red solid), McQuinn et al. (green dashed), 
FFRT (blue short-dashed), and FFRT-S (magenta dotted). Plot- 
ted is the likelihood that a given ionized patch is part of a region 
with size R. The ionized bubbles get noticeably larger with in- 
creasing ionization fraction with all the methods and the trends 
and peak scales agree among all the methods. 
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Fig. 3. — Ionization field power spectra for the four schemes. 
As reionization progresses and bubbles grow and merge, the power 
shifts to larger scales. The semi-numeric schemes have more power 
on large scales due to their larger and more connected bubbles. 
The small scale upturn in power at k > 8 h/Mpc is likely due to 
shot noise that is inherent in all of the schemes. 



advantage of being quickly computable with Fast Fourier 
Transforms. The bubble size distribution has a charac- 
teristic peak size that shifts toward larger scales as HII 
regions grow and reionization progresses. 

The distributions from the four schemes are generally 
in good agreement, as expected from the visual inspec- 
tion of the ionization maps. The two RT simulations 
show very good agreement on most scales and at all 
times. Some minor differences are found in the sizes of 
the largest and rarest of HII regions. At early times, 
the FFRT algorithm produces more small bubbles while 
the FFRT-S produces less relative to the RT simulations. 
These differences might be due to the disparate default 
filters. The top-hat filter of the FFRT-S scheme, oper- 
ating on relatively few discrete sources early in reioniza- 
tion, tends to produce less connected HII regions (and 
hence more small bubbles according to the definition 
used here). On the other hand, the sharp k-space filter 
used in FFRT can include ionizing photons from more 
distant sources through its non-vanishing real space fil- 
ter response, which should result in more connectivity 



of HII regions. Later in reionization when sources are 
more abundant, both semi-numeric schemes tend to pro- 
duce slightly larger bubbles because of more connectivity 
between ionized regions. 

5.2. Power spectra 

Another measure of reionization morphology is the 
auto power spectrum of ionized regions Pxx, which is 
relevant because it is directly related to the 21cm power 
spectrum (see Equation [2|). Figure |3] shows the evo- 
lution of t he ionization power spectra. As previ ously 
noted (e.g. iMcQuinn et al.ll2007btlZahn et al.ll2007[ ). the 
power shifts to larger scales as the bubbles grow and com- 
bine. Furthermore, the power spectrum is relatively flat 
as a function of scale during the final half of reionization, 
though semi-numeric studies predict a drop in power on 
scales larger than these simulation boxes (Mesinger & 
Furlanctto 2007). While all of our schemes show these 
trends, it is interesting to note that they are not as 
obvious ly present in the RT simulations of llliev et a"L1 
(|2006bD . 

The small differences in Pxx between the four schemes 
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Fig . 4. — Cross correlation coefficients between the Trac & Cen 
(2007) ionization fields and those from the other schemes. The 
two RT simulations show very good agreement on all scales at all 
times. Of the semi-numeric schemes, the FFRT -S comes clo s est to 
matching the RT results. For comparison, the Zahn et al] (120051 . 
cyan dot-dashed) scheme is shown only in the middle panel. It 
shows poorest correlation owing to usage of the un-evolved density 
field and poor tracing of the density by the spherical top-hat filter. 



are consistent with the differences seen previously in 
the bubble size distributions. At later times, the semi- 
numeric models have ~ 30% more large-scale power than 
the RT simulations because of their larger and more con- 
nected bubbles. As a more stringent test of the agree- 
ment between the ionization helds, we consider the cross 
correlation coefficient r Xl x 2 = PxtX 2 / \ / Px 1 x 1 Px 2 x 2 , 
where -PxiX 2 is cross power spectrum of ionization helds 
Xi and X2. While the schemes may agree on Pxx by 
having similar distributions of bubble sizes, rx x x 2 specif- 
ically tests whether the ionized cells are in the correct 
locations by comparing the phases of the Fourier modes 
of the transformed ionization helds. 

Figure |H s hows the cross correlation between the 
Trac fc Cenl (120071) io nization helds and those of 
McQuinn et all (|2007bf) . FFRT, and FFRT-S. On all 
scales, the ionization structure in the RT simulations 
agrees well, with rxiX 2 > 0.8 for all k out to the Nyquist 
mode. The slightly weaker correlation on scales near the 
grid spacing are due to small differences in the locations 
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Fig. 5. — Cross power spectra between the ionization and density 
fields for the four schemes. Sec Figure[6]for the corresponding cross 
correlation coefficients. 



and shapes of the ionization fronts. The semi-numeric 
schemes show progressively weaker correlation as reion- 
ization progresses, with rx ± x 2 > 0.8 at fc < 1 h/Mpc 
when Xi 0.5. On most scales, the FFRT scheme performs 
less well than the FFRT-S scheme, because the discrete 
nature of sources is not captured in this excursion-set 
based analytic scheme. Note the signihcant improve- 
ment in the cross correlation coefficient in the FFRT 
scheme compared to the original sem i -nume ric reioniza- 
tion scheme, discussed in iZahn et all (|2005l ) [and shown 
only in the middle panel]. This improvement can be at- 
tributed mostly to the usage of the true density field 
(rather than the linear field) to calculate the collapse 
fraction. 

We also investigate how the bubbles are distributed 
with respect to the large-scale overdensities of dark mat- 
ter and gas. As suggested already by many studies, large- 
scale overdense regions tend to be ionized earlier than 
large-scale underdense regions because they harbor more 
sources. The semi-numeric filtering techniques are moti- 
vated by this general description of the reionization pro- 
cess. In order to quantify the validity of this assumption, 
we calculate the cross power spectrum Pxd between the 
ionization and density fields, which is also directly related 
to the 21cm power spectrum (see Equation [2]). Figure [5] 
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Fig. 6. — Cross correlation coefficients between the ionization and 
density fields for the four schemes. The correlation is perfect on 
large scales, but it weakens with increasing k because the sources 
(and sinks) are not perfectly correlated with the density field on 
small scales. 



shows Pxd and Figu re \5\ shows the cross correlation coef- 
ficient txd = -Pxd / \J PxxPdd ■ In general, the ionization 
and density fields are perfectly correlated on large scales, 
but the correlation gets progressively weaker on smaller 
scales. This is expected as the sources (and sinks) are 
not perfectly correlated with the density field on small 
scales either. 

The two RT simulations sho w good agre e ment over- 
all in both P X d and r XD - The lTrac fc Cenl (|2007l ) sim- 
ulation appears to track the density field more closely, 
particularly for the xi = 0.25 case, with txd differing 
by ~ 0.1 on scales near the grid spacing. The FRRT- 
S scheme compares very w ell with the tw o RT s imu- 
lations, especially with the iMcQuinn et al.l (|2007b[ ) re- 
sults. There are some minor differences on large scales, 
but the good agreement overall validates the assump- 
tion discussed above. Not surprisingly, the FFRT scheme 
produces larger values of Pxd and txd (particularly at 
fc > 1 h/Mpc) compared to the other schemes because in 
this case, both the sources and ionization are determined 
by the evolved density field. Note that radiative transfer 
effects and poisson fluctuations in the source distribution 
tend to decrease the cross correlation. 



6. CONVERGENCE IN THE 21 CM SIGNAL 

One of the main upcoming applications of reionization 
models will be to aid in the interpretation of upcoming 
high-redshift 21cm observations. Hence, we briefly ex- 
plore the differences in the 21cm signal predicted by our 
various schemes. We only study the power spectrum as 
this is a simple, often-studied statistic, which can encode 
infor mation about the reionization state of the Universe 
(e.g- ILidz et~all[2008h 

Ignoring redshift-space distortions, the excess 21 cm 
brightness temperature over that of the cosmic mi- 
crow ave background from em ission at a redshift of z is 
(e.g., IZaldarriaea et al.l 120041 ): 



5T « 26(1 + 5)x H i 
0.15 



- CMB 



1 + z 



10 



T s ) V 0.022 

1/2 

mK 



where 5 is the density contrast of gas, Ts is the spin 
temperature, Pcmb is the temperature of the cosmic mi- 
crowave background, and xhi = 1 — Xi_v is the average 
neutral fraction. The first-order expansion in ionization 
and density fluctuations for the 21 cm power spectrum 
is 



Pnc m (fc)=X- 6 2 [Pxx-2x HI P : 



XD 



(2) 



where P& is the average brightness temperature in regions 
with xhi — 1, and P is the power spectrum of the ion- 
ization (X) and density (D) fields. Equation ^ ignores 
peculiar velocities and assumes Ts ^> Pcmb- 

In Fig. [3 we show the 21 cm power spectra 
computed using the four algorithms. We assume 
Ts ^ Pcmb, which previous studies have concluded 
should be accurate throughout most of reioniza tion (e.g., 
lFurlanettoir2006l : iPritchard fc Furlanettoll2007l) . We also 
neglect redshift-space distortions, which mainly affect 
smaller scales than the bubble scale during reioniz ation 
(|McQuinn et al.l 120061 : iMesinger fc Fur lanettol 120071 ). 

As anticipated from the agreement between predictions 
for Pxx and Pxd, the prediction for P2i C m between the 
two radiative transfer simulations is in good agreement, 
differing by < 10% on all scales. The FFRT-S tends to 
have more power by ~ 30 percent compared to the RT 
simulations, though this difference decreases as reioniza- 
tion progresses. The FFRT algorithm has a more com- 
plicated behavior, tending to over-predict the large-scale 
power, and under-predict the small-scale power. The 
over-prediction of power on large scales is likely due to 
its over-connected HII regions. The under-prediction of 
power on small scales is due mainly to its overprediction 
of Pxd (see Figure [5] and eq. [2|) . Nevertheless, the 21 
cm power spectra from all of the schemes agree with one 
another at the 10s of percent level. 

The error bars in the middle panel of Figure [7] are the 
forecasted precision that the MWA and the SKA can 
measure P^rm, assuming the configurations described in 
IMcQuinn et alJ (|2006l ). These errors assuming a 1000 hr 
integration and a bandwidth of 6 MHz. The differences 
between the four schemes are smaller than the MWA 
1-cr errors at fc > 0.5 h/Mpc. On larger scales, how- 
ever, uncertainty in the semi-numerical modeling might 
bias reionization constraints within the errors of MWA, 
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Fig. 7. — Comparison of the 21 cm power spectrum in the four 
schemes, ignoring peculiar velocities. At t he intermed iate ioniza- 
tion fraction we also compare the results of Zahn et al. ( 2005|) and 
Mesingcr & Furlanctto (2007), which agree less well with the sim- 
ulations than the improved schemes proposed in this paper. Also 
included are the MWA and SKA forecasted errors. 



and especially with the SKA. However, it is debatable if 
we can accurately predict either the 21cm power spec- 
trum or the relative differences among the ionization 
schemes on these large scales, given our limited box 
size. 18 Moreover, the relative differences among the ion- 
ization schemes are much smaller than the differences 
which can result from some astrophysical uncertainties 
during the EoR , such as the galaxy ioniz ing-photon lumi- 
nosity function iMcQuinn et al.l ()2007bD and the number 
and distribution of absorption systems (|Furlanetto fc Ohl 
[200l IChoudhurv et al.|[2009h . 

7. CONCLUSIONS 

We have compared the results of different algorithms 
for computing the morphology of reionization. Two 
of these algorithms were radiative transfer simulations, 
which made disparate approximations for propagating 

18 Since all of the power spectra are computed from one realiza- 
tion of the density field, it is unclear how the mean signal and the 
differences between the models would change if the power spectra 
were computed from a much larger simulation box, or from an en- 
semble average (using different realizations of the initial conditions 
sampling the modes larger than our box size) . 



HII fronts from the sources, and the other two are 
less computationally expensive semi-numeric schemes, 
which use excursion-set formalism to identify ionized 
regions. All of these schemes (or closely related vari- 
ants) have b een used in publ is hed s t udies of reion- 
ization (e.g.. IZahn et al.l 120051 120071: IMcQuinn et alj 



2007 
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iMesinger fc Furlanettol 120071: IMcQuinn et all 
2007at iTrac fc Cenl 120071: iMesinger fe Furlanettol 12008 
Trac et al.ll2008|) . These algorithms were all compared 
on equal footing, using the same realization of the large 
scale matter distribution on which to generate the ion- 
ization field. 

The two radiative transfer schemes are in excellent 
agreement with each other (with the cross-correlation co- 
efficient of the ionization fields given by rxx > 0.8 for 
fc < 10 h/Mpc) and in good agreement with the analytic 
schemes (rxx > 0.6 for fc < 1 h/Mpc). When used to 
predict the 21cm power spectra, one of the most impor- 
tant applications of reionization simulations in the near 
future, all ionization algorithms agree with one another 
at the 10s of percent level. The differences between these 
schemes are smaller than other uncertainties with mod- 
eling the EoR, such as the properties of the sources and 
sinks of ionizing photons. This agreement suggests that 
the different approximations involved in the ray tracing 
algorithms are sensible and that semi-numerical schemes 
provide a numerically inexpensive, yet fairly accurate re- 
alization of the reionization process. The speed and accu- 
racy of these schemes can aid in efficient parameter space 
studies and in the interpretation of future observations. 
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8. APPENDIX: FURTHER TESTS OF THE FFRT SCHEME 
8.1. Distributions of the Ionization Fraction 

Here we investigate how well the analytic schemes per- 
form at computing the level of partial ionization in sim- 
ulation cells. This is where we expect the most disagree- 
ment between the RT and FFRT, as these cells contain 
the edges of HII regions/ionization fronts, where detailed 
radiative transfer effects are the most important. 

We remind the reader that the FFRT scheme assigns 
partial ionization fractions based on the collapse frac- 
tion inside each cell. Therefore, it takes into account 
the ionizing photons from sub-grid sources, which were 
not luminous enough to ionize the entire cell. This is 
done in anticipation of the use of FF RT on very-larg e 
scales, with fairly large cells sizes (e.g., iMesingen 120091 ) . 
As the cell size increases, fluctuations in the source num- 
ber among neighboring cells decreases, and our approx- 
imation becomes increasingly accurate. However, as the 
cell size becomes small (M ce u ~ Af mm ), partially ionized 
cells increasingly correspond to the edges of ionization 
fronts driven by external sources, instead of the sub-grid 
sources chewing away at their local HI as assumed by the 
FFRT scheme. Here we wish to quantify these limits, 
checking how well the FFRT reproduces the distribution 
of partially ionized cells, given a certain cell size. 

In Fig. [5] we show the fraction of simulation cells with 
a given ionization fraction, at (z,xiy)= (8.49, 0.25), 
(7.56, 0.51), (7.11, 0.72), top to bottom. Red circles, 
green triangles, blue squares, and magent a pentagons 
corres p ond to the ionization a lgorithms of iTrac fe Cenl 
(|2lMhlMcQuinn et al.l (|2007aD . FFRT, and FFRT with- 
out Poisson fluctuations in the halo number (see , re- 
spectively. The left panel of Fig. [5] shows the ionization 
fraction distributions for the 256 3 grid used throughout 
this paper. The two RT schemes agree extremely well, 
signaling again that these schemes have converged. Note 
that on this grid scale, and assuming a relatively soft 
stellar ionization spectrum, the ionization field can in- 
deed be well represented with a binary field. The FFRT 
algorithm tends to under-prcdict the number of mostly- 
ionized cells. This under-prediction is less severe in the 
early stages of reionization, when the surface area of ion- 
ization fronts is smaller. 

Another interesting thing to note from the left panel 
is that our fiducial FFRT model which includes Poisson 
scatter in the sub-grid source number, has no partially 
ionized cells with ionization fractions less than ~ 0.1. 
This is because a single source with our minimum halo 
mass M m i n is luminous enough to ionize its host cell to 
a degree > 0.1. On the other hand, if one ignores that 
sources are discrete and instead uses the mean condi- 
tional collapse fraction (magenta pentagons), the distri- 
bution of partial ionizations agrees fairly well with the 
RT in this regime. This indicates that partial ionizations 
in the 256 3 grid with its Ax = 0.56 Mpc cells is domi- 
nated by unresolved ionization fronts driven by sources 
outside of each cell, as discussed above. Hence, these cells 
are too small to justify the assumption of partial ioniza- 
tions driven by sub-grid sources, and the good perfor- 
mance of the FFRT without Poisson scatter at Xi < 0.1 
is coincidental. Fortunately the ionization field on these 
scales is well approximated by a binary (fully neutral or 



fully ionized) field, and so correctly computing the partial 
ionization fraction is not very important for observables 
such as the 21cm power spectrum. 

As discussed above, we expect our partial ionization 
scheme to be most effective when cell sizes are large, 
M ce ii ^> Mmin. Unfortunately, we do not have RT fields 
for very large scales and cell sizes with which to check 
this assumption. Instead we approximate larger cell sizes 
by smoothing the 256 3 density field down to 128 3 , and 
using this lower-resolution grid as the root grid for the 
FFRT algorithm. We compare the resulting distribu- 
tions to the same RT fields (i.e., with the 256 3 root grid), 
but smoothed down to 128 3 . The resulting distributions 
are shown in the right panel of Fig. [5J These distribu- 
tions are not as smooth as the 256 3 ones, due to smaller 
number statistics. However, we can note that the FFRT 
scheme does in fact perform better on these Ax = 1.1 
Mpc cells. This may or may not be a physical effect, due 
to the increasing relevance of sub-grid sources in setting 
the partial ionization fraction of larger cells. We certainly 
expect agreement to improve with larger cell sizes. 

Another interesting point evident from the left panel 
of Fig. is that not including the discreteness of sources 
on this cell size results in a notable lack of fully neutral 
pixels. This is to be expected since the mean value of 
the collapse fraction (which is the quantity directly com- 
puted in PS formalism) is never zero for M ce \\ > M m ; n . 
Therefore, including the discrete nature of sources, such 
as with our Poisson scatter approach, is important in 
identifying pristine, fully-neutral regions on this scale. 

8.2. Impact of filter choice and cell flagging algorithm 

Here we would like to explore how the filter choice 
and choice of cell flagging algorithm affect the re- 
sulting ionization fields. Note that these are two of 
the main differences between ou r FFRT, FFRT- S algo- 
rithms, and those introduced in iZahn et all (|2005l) and 
iMesinger fc Furlanettol (|2007l ). To do so we show the 
Pxx,Pxs, and P2icm power spectra of these models, 
together with the simulation that used the scheme of 
ITrac fc Ceil (|2007| ) for one redshift in Figure H 

These older schemes do not include partial ionizations, 
though this has a minimal e ffect on thes e scale s as seen 
above. Aside from this, the IZahn et"aLl f|2005T ) scheme 
is different from the FFRT mainly in using a top-hat in- 
stead of a sharp k-space filter and operating on the linear 
density field, instead o f the e volved density field. The 
IMesinger fc Furlane tto l|2007[ ) scheme is different from 
the FFRT-S mainly in painting the entire spherical re- 
gion enclosed by the filter as ionized ("flagging entire 
sphere" ) , instead of just the central filter cell ( "flagging 
center cell"). 

All models do a decent job of predicting the ionization 
power spectra. However, it can be seen from the mid- 
dle panel that both of the older schemes under-predict 
the cross-correlation of the ionization and the density 
fields on moderate scales. As shown in eq. ((2J, this 
underprediction results in an increase in 21cm power, 
as can be seen in the bottom panel. The sharp k- 
space filtering (FFRT) reaches more small scale features 
in the density field, comparable to the RT algorithm, 
conversely tending to over-predict the ionization-density 
correlation P xc i on small scales, yielding a slight under- 
estimate of P2icm there. Note that our FFRT-S scheme 
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Fig. 8. — The fraction of simulation cells with a given ionization fraction, at (z, Xiv)= (8.49, 0.25), (7.56, 0.51), (7.11, 0.72), top to 
bottom. Red ci r cles, gr een triangles, blue squares, and magenta pentagons correspond to the ionization algorithms of Trac & Cen (2007), 
IMcQuinn et al.l 12007a! ) . FFRT, and FFRT without Poisson fluctuations in the halo number (see [|4j, respectively. The left panel was 
created from the 256 3 grid (cell size Ax = 0.56 Mpc) used throughout this paper, while the right panel was created from a 128 3 grid (cell 
size Ax = 1.1 Mpc). Note that the distributions corresponding to the RT simulations in the right panel were computed from just the 
(boxcar filter) smoothed 256 3 grid. However the FFRT ionization fields were generated directly from the low-resolution 128 3 root grid of 
the low-resolution density field (i.e., the 128 3 FFRT ionization fields are not just smoothed versions of the 256 3 fields). 



does very well in predicting this cross-correlation, sug- 
gesting that one has to use the center-flagging procedure 
and either use k-spacc filtering or use discrete sources 
not to get this underestimate. Interestingly, the shape 
of the power spectrum is ac curately predicted by the 
IMesinger &: Furlanettol (|2007T ) algorithm. Overall the 
new semi-numeric schemes introduced here fare best in 
comparison to the radiative transfer simulation. 

To summarize, if one is interested in just the ioniza- 
tion fields, all models do a good job. The "flagging the 
entire sphere" algorithm performs slightly better than 
"flagging center cell" in "by-eye" and cross-correlation 



comparisons with RT simulations. The situation is more 
complicated if one is interested in the 21cm fields, as 
cross-correlations with the density field and higher-order 
terms become relevant (Lid z et al.ll2009f ) . In this regime, 
"flagging the center cell" (instead of the entire sphere) 
seems to predict a better cross-correlation of the ion- 
ization and density fields, and is a simpler and faster 
algorithm. To get accurate 21cm power spectra, one 
should additionally either use a sharp k-space filter on 
the evolved density field (FFRT) or use discrete sources 
(FFRT-S). 
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